PROCEEDINGS 

°F SCIENCE 



(N 

O 

(N 

D 
00 



(N 

> 
O 

<N 



Random Matrix Models for Dirac Operators at finite 
Lattice Spacing 



Mario Kieburg 

Department of Physics and Astronomy, State University of New York at Stony Brook, NY 
11794-3800, USA 

E-mail: jnario . kieburg@stonybrook . edu 



Jacobus J. M. Verbaarschot 

Department of Physics and Astronomy, State University of New York at Stony Brook, NY 
11794-3800, USA 

E-mail: |jv@chi . physics . sunysb . edu 



Savvas Zafeiropoulos 

Department of Physics and Astronomy, State University of New York at Stony Brook, NY 
11794-3800, USA 
E-mail: 



szafeiroSic. sunysb . edu 



X: 



We study discretization effects of the Wilson and staggered Dirac operator with Nc>2 using chi- 
ral random matrix theory (chRMT). We obtain analytical results for the joint probability density 
of Wilson-chRMT in terms of a determinantal expression over complex pairs of eigenvalues, and 
real eigenvalues corresponding to eigenvectors of positive or negative chirality as well as for the 
eigenvalue densities. The expUcit dependence on the lattice spacing can be readily read off from 
our results which are compared to numerical simulations of Wilson-chRMT. For the staggered 
Dirac operator we have studied random matrices modeling the transition from non-degenerate 
eigenvalues at non-zero lattice spacing to degenerate ones in the continuum limit. 
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1. Introduction 

Chiral RMT is a powerful mathematical tool to calculate eigenvalue correlations in the infrared 
limit of quantum chromodynamics (QCD) and has been successfully applied to Dirac spectra for 
almost two decades [|l|]. In the low energy limit QCD exhibits universal behavior that agrees with 
chRMT with the same symmetries as the Dirac operator. In the continuum limit as well as for 
finite lattice spacing, a / 0, random matrix ensembles have been constructed that reproduce the 
eigenvalue correlations of the lattice QCD Dirac operator. 

Discretization effects of the Wilson Dirac operator were already studied with the help of chiral 
perturbation theory in the p-limit ^ ^]. Furthermore it was shown that the chiral Lagrangian 
agrees with Wilson chRMT in this limit [||, ^, ^. Random matrix theory enables us to get results 
which were not accessible before. We derive the joint probability density (jpd) of the non-Hermitian 
version for the random matrix ensemble proposed in Ref. ^ and present its eigenvalue densities 
in the microscopic limit which have to agree with the low lying eigenvalues of the Wilson Dirac 
operator. In our representation we are able to distinguish between the eigenvalue densities of the 
complex eigenvalues and of the real eigenvalues with eigenvectors of positive and negative chirality. 

Staggered fermions are widely used in lattice simulations mainly because of their low com- 
putational costs. We study a matrix model that reproduces the single trace terms in the chiral 
Lagrangian (reduced to staggered fermions in two dimensions) incorporating taste breaking effects 
at order [16|. It is the model proposed in Ref. |15] for staggered fermions in two dimensions 
which we expect to be analytically solvable. 

In Sec. ^, we consider the Wilson-chRMT [^] and study its spectral properties. A random 
matrix ensemble for the staggered fermions is proposed in Sec. ^ 

2. Wilson fermions 

Wilson-chRMT. The ensemble introduced in Ref. [||] is 
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For a = this gives the chiral Gaussian unitary ensemble (chGUE) that describes the infrared Dirac 
spectrum of continuum QCD The complex n x (n + v) matrices W and preserve chiral 
symmetry, whereas the Hermitian matrices A and B break the chiral symmetry and are identified 
with the Wilson term. In the microscopic limit (n — )• oo) the rescaled lattice spacing = na^/2, the 
rescaled eigenvalues z = 2nz as well as the index of the Dirac operator v are kept fix. Then, spectral 
correlations of Wilson-chRMT become universal and agree with infrared Wilson chiral perturbation 
theory (xPT) The volume of space-time V is identified with the matrix dimension n. 

Let a be the physical lattice spacing. Then, the relation to the rescaled quantity is a = ^/WsVa. 
The low energy constants We and W^ corresponding to the squares of traces in the action of the 
Goldstone bosons [Q, |6|] are suppressed in the large Nc limit [||] and will not be considered here. 
We are interested in the quenched case and a more general setting will be considered elsewhere. 

Without loss of generality, let v > 0. Due to the 75 = diag(l„, — l„+v)-Hermiticity of Dw 
((75^w) ' = JsDw), the eigenvalues are either real or come in complex conjugated pairs. There are 
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V generic real eigenvalues corresponding to the v zero modes at a = 0. Moreover, 2/ (0 < / < «) 
additional real modes may appear when I complex conjugated eigenvalue pairs enter the real axis. 

Often the Hermitian version of the Wilson Dirac operator D5 = 75DW is studied because it 
simplifies lattice simulations [10] as well as RMT calculations ^ ^ However, only Dw 
is directly related to the chiral symmetry breaking which is our main motivation for studying the 
non-Hermitian version. 



The joint probability distribution. Analogously to Hermitian matrices, the 75 Hermiticity al- 
lows us to quasi-diagonalize the matrix Dw by a non-compact unitary matrix U G \J{n,n + v). 



Dw = and X^'^ = diag 





' X2 yi] 








y -yi X2 J 



(2.2) 



The real matrices x\,X2, yi and X3 are diagonal with dimensions l,n — l,n — l and Z + V, respectively. 
The ensemble decomposes in « + 1 disjoint sets differing in the number of complex conjugate 
pairs, n — l, or equivalently, in the number of real modes, 21 + v. The / complex pairs are given by 
{z2,z*2) = {x2 + ly 2,^2- ly 2)- 

Let Z = (zir, . . . ,z,jr,Zii, • • • ,Zn+v,i) G C^"+^ be the 2n + v eigenvalues of Dw where we ignore 
which of them are complex or real. The jpd is given by a sum over non-compact coset integrals 
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with the Haar measure J/Xg, on the coset G/ = \J{n,n + v) /[U"+'+^(l) x 0""'(1, 1)] and the Van- 
dermonde determinant A2„+v(2)- After the integration we are able to perform the sum by introduc- 
ing Dirac delta functions and find |jT3|] 
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The function erfc is the complementary error function and b^^\x^iy') = b{x)b{y). The 75- 
Hermiticity enforces a breaking of the permutation group 5(2?! + v) to S(n) x5(« + v) which 
is reflected in the expression ( |2.4| ). An expansion in the Dirac delta functions yields all « + 1 terms 



in the sum (2.3) corresponding to a fixed number of real modes. The two-point weight g2 splits 
into a weight for the real modes, g2x, and for the complex conjugated pairs, g2z- Such a structure is 



already known from the real Ginibre ensemble [11 ] and its chiral counterpart [12]. 
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Figure 1: These 3d-plots show for small lattice spacing (a = 0.1, left figure) and for a large one (a~ 1, 
right figure). The index is v = 1. Oscillations are well distinguishable for small a whereas they disappear 
for larger lattice spacing. Notice that a in the right figure is one order larger than in the left one. 



The eigenvalue densities. The integration over all eigenvalues except one yields the three eigen- 
value densities of real modes, Pr for positive chirality {{yIYsIy) > 0), pi for negative chirality 
((wIYsIy) < 0), and of complex pairs, Pc, i.e. 

I p{Z) n d[zj] = pAxir)5{y,r) + / P{Z) n d[zj] = Pi(A-ii)5(yn) + (2.6) 



In both integrals the determinant ( |2.4[ ) can be expanded either in the first row or in the first column. 
The densities Pc and p^ can be expressed in terms of the Nf = 2 partition function 

Pc{z) = g2c{z){z- z*)Zj^^=2{z.,z*;a), p,{x) = J dyg2rix,y){x-y)Z^^^2ix,y;a). (2.7) 

The eigenvalue distribution of the left handed modes contains an additional contribution from the v 
generic real modes which is the chiral distribution Px= Pi — Pr- To express p^, into known partition 
functions additional rows of some of the determinants have to be expanded. For v = 1 and v = 2 
we checked that this result agrees with previously derived expressions [Q, ^]. 



In the microscopic limit (n — )• oo) we obtain 013| ] 
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(2.8) 

(2.9) 
■d[s].{2.\0) 



The functions sine, erf, 7/, Ki and 5^'^ are the sinus cardinalis, the generalized incomplete error 
function (erf(Z?,c) = erf(c) — erf(Z7)), Bessel function of the first kind, the modified one of the 
second kind and the /-th derivative of the 5 function, respectively. We employed the abbreviations 
D(pk = sin^(((pi — (p2)/2)d(pid(p2 and Aj = 2a[cos(pj —x/Scf). Notice that only the singular part 
of Ki contributes due to the 5 function and that p^ vanishes for v = 0. 
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Figure 2: Left: Log-log-plot of the average number of additional real modes over a. The Monte Carlo 
simulations (symbols) confirm our analytical results (solid curves). The matrix dimension and the number 
of matrices vary in the simulations. Right: The distribution Preai is governed by at small a (green, 
a~l /-\/200) but is increasingly dominated by Pr (red, a—l/ \/T). We compare our analytical results (solid 
curves) with Monte Carlo simulations (histograms, 200000 matrices, n ~ 50, V = 1, bin size is 0.4). 

For sufficiently small a, Pc is broadened by a Gaussian along the imaginary axis, the oscilla- 
tions of the continuum limit are distinct and near the real axis it behaves like 3;^+^ When increasing 
the lattice spacing the oscillations disappear and the behavior near the real axis becomes indepen- 
dent of V and y. The distribution develops a box-like shape along the imaginary axis with 
support \x\ < 8S^. In Fig. ^ we show 3D-plots of pc. 

An important quantity to measure the effect of a finite lattice spacing is the average number of 
the additional real modes. It is given by [ pl| ] 



add 



2 / p. 



dx - 



2k 



-Aa sin 
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a, a » 1. 



(2.11) 



At small lattice spacing contributions for non-zero index are suppressed whereas Nadd for large a is 
independent of the index. This is shown in Fig. ^. 

At small lattice spacing, the distribution of the real eigenvalues, p^ai = Pr + Pi = 2pr + P;^, is 
approximated by the v-dimensional Gaussian unitary ensemble. The width is 2a in this regime. 
For increasing lattice spacing the support increases to \x\ < 8S^ (see Fig. |2|) and develops a square 
root singularity at the edges, i.e. Preai(-''^) = l/[(27r)^/^a| + v/[K^y {Scf-y —x^]. The first term is 
due to Pr while the singularity comes from p^. 



3. Staggered fermions 

In this section we introduce an ensemble to study the taste breaking of the staggered Dirac 
operator. For simplicity we only consider the case of two tastes with the ensemble given by 



£>st 



W + aiW a2C 



P{D t) g-«(tr'*V2+trBV2+trwwt+trw/ivt) 



(3.1) 

where W,W aienx (n + v) and B, C are (« + v) x + v) and {n x n) complex matrices respectively. 
The matrix dimension n is identified with the spacetime volume V as in the case of the Wilson Dirac 
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Figure 3: Left: The spectral flow of staggered eigenvalues for v = 2 as a function of a using a fixed 
realization of the random matrix ensemble (3.1) for n ~ 50. The first avoided level crossing determines the 
perturbative regime. Right: The quenched spectral densities of Dst for values of a (nonrescaled) and v given 
in the legend of the figure obtained by Monte Carlo simulations with « = 50 and 100000 matrices. In both 
figures we set fli ^ a2 — a. 



operator. For a = this model has two flavors with v zero modes for each flavor while at a 7^ 
with interacting tastes, the zero modes are absent and the continuum SU(2) flavor symmetry is 
broken. A more general random matrix model with four tastes and additional taste breaking terms 



was introduced in Ref. [15] 



Studying the ensemble (3. 1 ) in the limit, « — )• 00, at rescaled quark masses M = InM and lattice 
spacing = allows us to obtain universal results for the eigenvalue correlations which will 

be worked out elsewhere. Here, we evaluate the partition function of A^t fermions corresponding to 



the matrix model ( |3.1[ ) which is given by the unitary matrix integral 

Zoc /" t/[/det^[/exp [a?tr[T3[/T3[/'^]+5-2tr[(Tif/)^ + (T2f^)^ + cc.] -trM([/ + [/t)l . (3.2) 

J\](Nf) L J 

It agrees with the e limit of the staggered chiral Lagrangian corresponding to the taste breaking 



pattern of the ensemble (p. l]), and is a special case of the result derived in Ref. []15[]. 

The spectral flow as a function of a shows avoided level crossings when lattice artifacts start 
dominating the Dirac spectrum (see Fig. || (left)). For small a we observe a linear behavior that 
follows from perturbation theory. At a = all eigenvalues have a degeneracy of two. The matrix 
W in Eq. ( p?l| ) lifts this degeneracy but does not give repulsion between the two subspectra. On the 
other hand, the matrices B and C lift the degeneracy and cause a repulsion between the subspectra. 

For cij <^ 1 the spectral density exhibits a chGUE with twice the number of flavors as for 
ay » 1. At a2 = the distribution of the zero modes is given by a Dirac delta function with weight 
2v at zero. For increasing a2 the distribution of these modes broadens gradually. Because they 
repel each other as well as the other modes, we always find a vanishing eigenvalue density at zero. 

4. Conclusions 

For the Wilson-Dirac operator the discretization effects become strong when a ^> 0.5, i.e. a » 
1 / (2^/WsV). Then characteristics of the continuum limit like the oscillations in the spectral density 
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or its behavior near the real axis disappear. In the limit of large a, the eigenvalue distributions 
become independent of the index v, have a support in \x\ < SWgVa^ and the distribution for the 
complex as well as for the eigenvalues of the right handed modes develop plateaus. The chirality 
distribution for large lattice spacing has square root singularities at the edges. 

At sufficiently small lattice spacings the low energy constant Wg can be extracted by combining 
the width of the Gaussian which broadens pc in the x direction, a = 2ay/W^/V /£, and the spacing 
of the projected eigenvalues onto the imaginary axis, AA = 7i/LV, yielding a /AX = 2a /n. Fur- 
thermore, additional real modes are highly suppressed by non-zero index v and, thus, are not much 
of a problem in lattice QCD since most configurations have an index |v| > for large volumes. 



The matrix model proposed above is the version with two tastes of the one in Ref. [15] and 
describes the main features of the low-lying spectrum of the two dimensional staggered Dirac 
operator. The scale 1 / y/n (i.e. 1 / ^/V ) determines the appearance of lattice artifacts with increasing 
lattice spacing a. We expect that it is analytically solvable like the ensemble for the Wilson Dirac 
operator. 
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